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Abstract 

The semi-Lagrangian method using the hybrid-cubic-rational interpolation function 
[M. Ida, Comput. Fluid Dyn. J. 10 (2001) 159] is modified to a conservative method 
by incorporating the concept discussed in [R. Tanaka et al., Comput. Phys. Com- 
mun. 126 (2000) 232]. In the method due to Tanaka et al., not only a physical 
quantity but also its integrated quantity within a computational cell are used as 
dependent variables, and the mass conservation is completely achieved by giving 
a constraint to a forth-order polynomial used as an interpolation function. In the 
present method, a hybrid-cubic-rational function whose optimal mixing ratio was 
determined theoretically is employed for the interpolation, and its derivative is used 
for updating the physical quantity. The numerical oscillation appearing in results by 
the method due to Tanaka et al. is sufficiently eliminated by the use of the hybrid 
function. 

Key words: Numerical method. Conservative method, Semi-Lagrangian, 
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1 Introduction 



The semi-Lagrangian method is a category of numerical methods for par- 
tial differential equations including an advection (or convection) term, and 
has been used mainly for atmospheric and geographic problems [1,2]. In the 
semi-Lagrangian approach, the advection of a physical quantity is solved as a 
transportation and interpolation problem on fixed grids. The trajectory of the 
Lagrangian invariant that will reach a grid point is computed backwardly by 
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the time integration of the advcction velocity, and the invariant at the depar- 
ture point, which is determined by an interpolation within a computational 
cell, is transported to the grid point. One can use a CFL-free time step in this 
approach. Since there are some well-established methods for calculating the 
trajectory such as the Runge-Kutta methods, recent studies on this approach 
are focused on the construction of interpolation functions. 

One of problems of this approach is the difficulty of constructing a completely 
conservative scheme. There are only a few studies attacking to this problem. 
In Ref. [3] Priestley proposed a quasi-conservative semi-Lagrangian method by 
coupling a flux- corrected transport method and a high-order interpolation. In 
Ref. [4], Leonard et al. proposed a conservative, explicit algorithm achieving 
stable high-CFL computations, by using an integral variable of a Lagrangian 
invariant as the dependent variable, and showed results for one-dimensional 
advection at constant velocity by coupling the algorithm with many interpola- 
tion schemes. In Ref. [5], Manson and Wallis modified the QUICKEST scheme 
[6] so that the stable high-CFL computation is achieved. Applications of their 
method to pure advection problems in non-uniform flow flelds are shown in 
Ref. [7]. In Ref. [8] Lin and Rood discussed the use of readymade Eulerian 
schemes (the MUSCL [9] and PPM [10] schemes) modified to achieve stable 
computations even with a large time step. 

Recently, Tanaka et al. proposed a way to solve this difficulty. In Ref. [11,12], 
they modified the CIP method, which is a non-conservative semi-Lagrangian 
method using the Hermite cubic interpolation function [13,14], into a conser- 
vative one by employing the cell-integrated quantity of a Lagrangian invariant 
as an additional dependent variable. A 4th- or 2nd-order polynomial is used as 
an interpolation function, and is constructed using a physical quantity, its first 
spatial derivative (only for the 4th-order one), and its cell-integrated quan- 
tity. The conservation of total mass is achieved by updating the cell-integral 
quantity in a conservation sense. The CIP method and its conservative variant 
using the 2nd-order polynomial (the CIP-CSL2 method) are briefly reviewed 
in the next section. 

The aim of this paper is to propose a variant of the CIP-CSL2 method by 
employing the hybrid-cubic-rational (HCR) interpolation function introduced 
recently by the author [15]. The hybrid function is constructed with a com- 
bination of the cubic polynomial and a rational function, and is proposed to 
obtain oscillation-free results having higher resolution than those given by 
the rational method [16] being a variant of the CIP. In the results using the 
CIP-CSL2 method one can sometimes observe numerical oscillation caused 
by the use of the classical polynomial [12]. As is well known, the numeri- 
cal oscillation becomes a serious problem for a certain examples such as the 
propagation of shock waves and the multiphase flow problems including large 
density jump. In Ref. [17] the use of the rational function employed in the ra- 
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tional method is suggested to overcome the numerical oscillation. The author 
expects that combining the hybrid function with the concept of the CIP-CSL 
can derive an oscillation- free, conservative semi-Lagrangian method that has 
higher resolution than that by the rational method. In Sec. 3 we show that the 
hybrid method can easily be modified into a conservative one, and in Sec. 4 we 
demonstrate the accuracy and stability of the present method. In this paper, 
we discuss only the case of CFL number < 1; however, we believe that the 
present method is applicable to a high-CFL condition by incorporating the 
concept discussed in Ref. [12]. 



2 CIP method and its conservative variants 



The CIP method [13,14] is a semi-Lagrangian solver for an advection equation. 



9f ^ n 



where / is a dependent variable and u is the advection velocity. The fact that 
the solution of this equation can be represented as 

f{x,t + At)^f{X{x,t),t) (2) 



allows us to solve Eq. (1) as an interpolation problem. Where X is the trajec- 
tory of fluid particle that locates at x at the time t + At, 

t-At 

X{x,t)=x+ J u{X{x,T),T)dT. (3) 



In the CIP method, the dependent variable / is interpolated by the Hermite 
cubic function, 

C{X) = f, + d,{X - X,) + C2^{X - x,f + CUX - x^)\ 



constructed using both and its derivative di ( = df/dx\i) defined at each 
grid points Xi, where i = 1,2, ■ ■ ■ , N, N is the number of grids, and C2j and 
C3j are determined from the condition of continuity [13]. Equation (1) and its 
spatial derivative, i.e., 

df{x,t + At) ^ dXix,t) df{X{x,t),t) 

dx dx dX ' ^ ' 
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arc used to update fi and di, respectively. In the case where we assume the 
first-order accuracy in time for example, Eqs. (2) and (4), respectively, are 
reduced to 



f{x,t + At) ^ f{x -u{x,t)At,t), 



df{x,t + At) 
dx 



1 - ^}<^At 
dx 



df{x-u{x,t)At,t) 
dX 



(5) 
(6) 



where At is the time interval during a computational step. In the case of 
Ui < 0, those equations are represented using the Hermite cubic function as 



(7) 

(8) 



where 



C2r = -(C+i + 2C-3^r+i/2)A, 

car = «+i + c - 25r+i/2)//^', 



on 



(9) 
(10) 



h is the grid width, and the superscript n denotes the time t — nAt. 
Basically, the CIP-CSL2 method is a solver for the equation of continuity. 



In this method, a quadratic function of 

Qi{X) = fi + 2 qli{X -Xi) + 3 q2i{X - Xif 



[12) 



is used to interpolate the physical quantity /, where the coefficients qli and 
q2i are determined, for example in the case of < 0, from the constrains of 



Xi + l 



Qi{xi+i) = fi+i and j Qi{x) dx/h = Pi+i/2, 
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as 

gl, = -(/,+! + 2/, -3a+i/2)A, (13) 
?2, = (/,+l + /,-2p,+l/2)//^^ (14) 

where p is the cell-integrated average of / within a computational cell, defined 
at the cell centers, and is used as a dependent variable in this method (We 
should note that the definition of p is shghtly different from that in Rcf. [11]). 
p is updated so that the total mass is conserved. The mass fiux flowing out 
from the cell [xj, Xj+i] through the point Xi during At can be calculated by 

Apih= I Qi{x) dx = M + qU'' + q2i^\ 



where we assume again the first-order accuracy in time. Using this quantity, 
p is updated as follows: 

PlV\/2 = Pm/2 + Apr+i - Apr- (15) 



The use of this formula allows us the complete conservation of the total mass, 
i.e., the total sum of p. 

The quantity / is updated by a time splitting technique [11]. We solve Eq. (11) 
by splitting into two phases as 

dt dx 
dt dx 



The former one shows the advection of / and, therefore, is solved as an inter- 
polation problem by using Eq. (12) as 



The latter one represents the change of / due to the compression or expansion 
and is solved in general by a finite difference technique [11,13]. 

As was pointed out already [12], the coefficients of the first- and second-order 
terms in Eq. (16) correspond to those in Eq. (8) by replacing / and p in Eq. (16) 
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with d and S, respectively. Furthermore, Eq. (15) also can be rewritten into 
a form like the conventional CIP. By introducing a variable D defined as 




one obtains 



(17) 



and 



PlVi/2 = {D\ 



in 



D-)IK 



where we assume Apg = 0, namely no transfer of mass occur at the bound- 
ary xq. Those equations reproduce Eq. (15), and replacements of D and / 
in Eq. (17) with / and d, respectively, yield Eq. (7) (Note that the former 
replacement reduces p to 5"). This result suggests us that the variants of the 
CIP method such as the rational method [16] and the HCR method [15] can 
be modified to a conservative one only by those replacements. 



3 A conservative method by the hybrid-cubic-rational interpola- 



In Ref. [15], the author proposed a numerical solver for the advection equation 
by employing both the cubic [13] and rational [16] functions. In the method, 
the following combination of those functions is used as an interpolation func- 
tion: 



where R and C show the rational and the cubic functions, respectively, and 
a denotes the weighting parameter whose range is a e [0, 1]. Equation (18) is 
reduced to the rational function for a = 1 and to the cubic one for a = 0. The 
value of a is determined theoretically so as to be the minimum of the values to 
which the convexity-preserving condition [16] is satisfied. Where the convexity- 
preserving condition is expressed as that Fxx{X) > for a concave data of 
di < S'j+i/2 < di+i or Fxx{X) < for a convex data of di > 5'j+i/2 > is 
satisfied for the region of X e [a;j,a;j+i]. The resulting formulae are 



tion 



F{X) ^aR{X) + {l-a) C{X), 



(18) 



F{k)^fi + dM + {Gli + G2,)e, 



(19) 
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^ ^d, + [Gl.^i^ + + (1 - am - A)]^, (20) 
ox Dj n 



where 



Di^Qi + {Pi - Qi)k, 
M,(M,-2) 



a 



Mi{Mi-2) + V 



Mi = max[2, max(^, ^)] 



with 



Pi — {Si+1/2 — di)h, Qi — (rfj+i — Si^i/2)h, 



and 



k = ^/h= -UiAt/h. 



For the case of -Uj > 0, we need replacements ofi + l^i — 1 and h ^ —h. 
This interpolation function provides oscillation-free results, unlike the cubic 
one, and higher resolution of solution than that by the conventional rational 
function [16]. Same as the conventional CIP, this method is constructed using 
a physical quantity fi and its first spatial derivative df, thus this method may 
be modified into a conservative method only by the replacements as discussed 
in the last section. Though the explicit definition and use of the cell-integrated 
average p may be more convenient in an actual application, we now use D, 
instead of p, and adapt the replacements of 



(21) 



for demonstration. Namely, the subroutine "HCR(/, (i)" of the hybrid method 
whose parameters are / and d is used as a function of D and /, i.e., HCR(D, /). 
The initial condition is set as follows. At first, / is initialized using an arbitrary 
function G{x), which expresses the initial profile of /, as 

f--G(xi). 
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Next, the initial value of D is calculated using a recurrence procedure of 

= 0, 

A° = DU + h{f^ + /°_i)/2, for z = 1, 2, • • • , TV. (22) 



Here the last term of Eq. (22) means the approximated integration of / over 
a cell between Xi-i and Xi. 

In the next section we show some numerical results using those procedures. 



4 Numerical experiments 



Example 1: 



The first example is the linear propagation of square waves. In this example 
we use u — 1, h — 1 and 



G{xi) = 



-1, for 13 < i < 21, 
1, for 40 < i < 48, 
0, elsewhere. 



Thus, the width of the square waves is 9/i. Figures 1 and 2 show the results 
using CFL = 0.2 at n = 200 and n — 2000, respectively. Here we show results 
using four types of method, the HCR, CIP (the case of a = 0), rational 
(a = 1), and modified rational methods. In the modified rational method [18], 
a switching technique, which is represented as 



7= \ 



1, di ■ di+i < 0, 
0, otherwise. 



is added to the conventional rational method to select the interpolation func- 
tion. The rational function is adopted in the case of 7 = 1 and the cubic 
one in the case of 7 = 0. It was pointed out that this technique sometimes 
breaks the preservation of the convexity of solution while this improves the 
dissipation property of the rational method [15]. In the present example, this 
technique is adapted by using /, not by rf, because of the replacements. In 
the figures, we see clearly that the HCR method provides most accurate re- 
sults, which are less diffusive than those by the conventional rational method 
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CFL = 0.2; STEP = 200 




10 20 30 40 50 60 10 20 30 40 50 




10 20 30 40 50 60 10 20 30 40 50 

X X 



Fig. 1. Linear propagation of square waves by (a): the HCR, (b): CIP, (c): rational, 
and (d): modified rational methods with the replacements of Eq. (21). Those results 
are with CFL = 0.2 at n = 200. The circles and the solid lines show the numerical 
and theoretical results, respectively. 



CFL = 0.2; STEP = 2000 




10 20 30 40 50 60 10 20 30 40 50 




10 20 30 40 50 60 10 20 30 40 50 60 

X X 



Fig. 2. Same as Fig. 1 but n = 2000. 

and less oscillatory than those by the CIP and modified rational methods. In 
Table 1, we show the calculated values at n = 200 around the point i = 12 
at which the left discontinuity of the negative pulse locates. The calculated 
data by the CIP and modified rational methods are positive near the left side 
of the discontinuity and become minimums at i = 14, which are smaller than 
— 1. On the contrary, the data by the HCR and conventional rational meth- 
ods decrease monotonously as i increases, and are not less than —1. Those 
results mean that only the HCR and conventional rational methods can keep 
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CFL = 0.2; STEP = 2000 




10 20 30 40 50 60 10 20 30 40 50 
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-0.5 
-1 
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10 20 30 40 50 60 10 20 30 40 50 60 

X X 



Fig. 3. Same as Fig. 2 but the replacements of Eq. (21) are not adapted; namely 
those schemes are solved as non-conservative methods. 



Tab: 



i 


Present 


CIP 


Rational 


Mod. Rat. 


4 





-0.000014 








5 





-0.000986 


-0.000004 





6 





-0.001887 


-0.000037 


0.000019 


7 


-0.000001 


0.004413 


-0.000345 


0.000223 


8 


-0.000044 


0.024674 


-0.002738 


0.000387 


9 


-0.001716 


0.032729 


-0.018069 


-0.000109 


10 


-0.052075 


-0.052964 


-0.09232 


-0.04052 


11 


-0.305191 


-0.304522 


-0.32068 


-0.277075 


12 


-0.681895 


-0.665011 


-0.67146 


-0.647653 


13 


-0.954887 


-0.955764 


-0.90682 


-0.951461 


14 


-0.999656 


-1.058063 


-0.982846 


-1.059713 


15 


-0.999996 


-1.029841 


-0.997188 


-1.031324 


16 


-0.999997 


-0.999959 


-0.998634 


-1.000665 


le 1 



Calculated values of the results corresponding to those shown in Fig. 1. 
the convexity of solution. 



For comparison, we shall show results without the replacement of Eq. (21). 
Figure 3 shows the results at n = 2000 given by using the same initial condition 
as the previous one. In this case, we set the initial value of the first derivative 
as d° = at all points. As was pointed out in Ref. [12], the CIP-CSL2 and CIP 
methods (the CIP with and without the replacements, respectively) give quite 
similar results each other. However, another methods provide obviously differ- 
ent results by whether the replacements are employed or not. It is interesting 
to point out that the result by the HCR method with the replacements is less 
diffusive than that by without the replacements, while we cannot explain this 
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result for the moment. 



Example 2: 



Next, we solve the problem used in Ref. [16] to compare the performance 
of the conventional rational method with those of the MUSCL [9] and PPM 
[10] schemes (See Fig. 1 in the paper). In this example, we use h — 1 and 
CFL = 0.2. The initial value of / is set to 



{i - 20)/ll, for 20 < i < 31, 
l-(i-31)/20, for 31 < i < 41, 
1/2, for 41 < i < 60, 

1, for 60 < i < 80, 

0, elsewhere. 



and the velocity, assumed as constant in time and space, is set to m = 1. 
Figures 4 and 5 show the results at n = 440, given by the four types of 
methods, with and without the replacements, respectively. In Ref. [16] it was 
concluded that the conventional rational method gives better representation of 
the sharp corner at the top of the triangular wave than those by the MUSCL 
and PPM methods. From the present result, however, we know that the HCR 
method and its conservative version achieve still better representation of that 
than the conventional rational method. The calculated maximum values of / 
at the corner are 0.935 (conservative HCR), 0.937 (HCR), 0.916 (conservative 
rational), and 0.923 (rational). Furthermore, we know that the conservative 
HCR method gives almost equivalent (or slightly better) resolution of the 
square wave to that by the PPM, which gives the best one among those by 
the schemes discussed in Ref. [16]. 



Example 3: 



Now we show results to a nonlinear problem in which the inviscid Burgers' 
equation, 

du du ^ , „, 



is selected as the governing equation [11,12]. This equation can be solved di- 
rectly by the conservative methods discussed in this paper. For the calculation 
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Fig. 4. Linear propagation of triangular and square waves. Those results are with 
CFL = 0.2 at n = 440. The circles and the solid lines show the numerical and 
theoretical results, respectively. 



CFL - 0.2; STEP - 440 



/ 




20 40 60 80 100 20 40 60 80 100 




20 40 60 80 100 20 40 60 80 100 

X X 



Fig. 5. Same as Fig. 4 but the replacements are not adapted. 

of we can use the procedures for / because Eq. (23) is an advection equa- 
tion (here / means that after the replacements). For the calculation of D 
{— J udx), we need to rewrite Eq. (23) into a conservation form as 



From this formula we found that the tranportation velocity u used in the 
procedures for D should be replaced with u/2 [11]. 
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Fig. 6. Result of inviscid Burgers' equation at t = 100 by (a): the HCR method with 
the replacements, (b): the CIP-CSL2 method, and (c): the HCR method without 
the replacements. The parameters are h = 1 and At = 0.1. The solid lines show 
results by the first-order upwind method with h = 1/5 and At = 0.1/5. 




Fig. 7. Maximum and minimum values in the results of inviscid Burgers' equation, 

as a function of time. Those results are by (a): the present and (b): CIP-CSL2 
methods with h = 1 and At = 0.1, and by (c): the first-order upwind method with 
h = 1/5 and At = 0.1/5. 
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The initial condition and the grid width are set to 

G{x) = 0.5 + 0.4cos(27ra;/100) 

and h = 1. Figure 6 shows results at t = 100 with At = 0.1. In this figure, 
the solid line shows the result by the first-order upwind method using the 
conservation formula (24) with h — 1/5 and At — 0.1/5. Same as the re- 
sult by the CIP-CSL2, that by the conservative HCR method represents the 
correct position of the shock wave. Figure 6(c) shows the result by the HCR 
method without the replacements. This result is obviously incorrect. Those 
results prove that the HCR method has gained complete conservation by the 
replacements. Figure 7 shows the maximum and minimum values of it as a 
function of time. In the result by the CIP-CSL2 wc can observe unphysical 
strong oscillation due to the numerical dispersion of a method using a classi- 
cal polynomial interpolation; however, such the oscillation cannot be seen in 
the result by the present method. This result allows us to consider that the 
present method is applicable to a shock problem even without any artificial 
viscosity, and is oscillation free even for the nonlinear problem. 

Example 4-' 

Now we shall try an experimental study which uses further replacements. 
We already know that the CIP method provides an almost equivalent result 
whether the replacements (21) are adapted or not. If the variables D and / 
are replaced with those integrated once more, what result is obtained? Figure 
8 shows results to Example 2 obtained by replacing and with 



and by updating those variables with the HCR and conventional methods. The 
parameters are same as those used in Example 2. fi shown in Fig. 8 is given 
by fi — (-E'i+i — S-Bj -I- Ei_i)/h^. Despite the adaptation of two times of the 
replacements, the CIP method provides an almost equivalent result to those 
given without or with one-time replacements. However, one can observe some 
differences in the results using another methods from the previous ones. The 
numerical oscillation around the discontinuities of the square wave appears 
even in the results obtained by using the HCR and the conventional rational 
methods. In the result by the modified rational method, the oscillation can be 
observed not only at the top of the square wave but also at the bottom of it. 




lO 
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CFL= 0.2; STEP = 440 





(b) CTP 
(2linie repl.) 



1 

0.8 








0.6 








0.4 








0.2 


/ (d) Mod. Rat. 
















40 60 

X 



Fig. 8. Same as Figs. 4 and 5 but the two-time replacements are adapted. 
5 An additional comment 



It is known that when no artificial viscosity is adapted, a numerical method 
for an advection equation based on the non-conservative formula cannot give 
correct propagation speed of shock in the region where the direction of the 
advection velocity changes spatially. In Ref. [11] the use of high-order poly- 
nomial that covers two cells are discussed for solving this problem. Now we 
present an alternative method to overcome this difficulty. 

Here we use again the Burgers' equation (23) as the governing equation. In- 
troducing U = u + C , where C is a constant, reduces Eq. (23) to 



dU_ jjdU_^dU 
dt dx dx 



0. 



This is also the Burgers' equation, but has a linear advection term (the last 
term of left hand side). We can solve this equation by splitting into two stages. 
First, we solve dU/dt + UdU/dx = and update U during the time interval 
At. Next, we shift the position of U by —CAt. If C is determined so that 
C/j > or C/j < for all i, we may be able to calculate the position of shock 
correctly. In the case where we set C so that the condition of m |CAt| = h 
(where m is a positive integer) is satisfied, the shift step can be solved only 
by one substitution like Ui-^i = Ui per m computational steps. This procedure 
is convenient since we need no additional interpolation function or difference 
equation; however, further discussions should be required for the practical 
application. 
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6 Conclusion 



In this paper we have derived a conservative variant of the hybrid-cubic- 
rational method for an advection equation. It was proved by some numerical 
experiments that the hybrid method is oscillation free even when it is modified 
to a conservative one. Interestingly, the conservative hybrid method gives less- 
diffusive representation of the discontinuities of the square wave than that by 
the non-conservative one. We should try to make this result clear theoretically. 

Wc expect that combining the present method with the hybrid interpolation- 
extrapolation method, which realizes the discontinuous representation of the 
density interface [19], yields a powerful tool for a semi-Lagrangian computation 
of multiphase flows like those treated in Refs. [20,21]. 
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